deficit_Model_Plots
Amani Chehimi
2025-05-12
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(knitr)
library(ggplot2)
library(dplyr)
library(lubridate)
library(MMWRweek)
library(gridExtra)##
## Attaching package: 'gridExtra'
##
## The following object is masked from 'package:dplyr':
##
## combinelibrary(readxl)
library(excessmort)
library(plotly)##
## Attaching package: 'plotly'
##
## The following object is masked from 'package:ggplot2':
##
## last_plot
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following object is masked from 'package:graphics':
##
## layoutlibrary(patchwork)
library(scales)##
## Attaching package: 'scales'
##
## The following object is masked from 'package:purrr':
##
## discard
##
## The following object is masked from 'package:readr':
##
## col_factorlibrary(cowplot)##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:patchwork':
##
## align_plots
##
## The following object is masked from 'package:lubridate':
##
## stampsource("~/Desktop/Thesis Work/Hospitalization/11-26 data/Deficit_Model_1.0.R")allcause_state_byweek <- read_excel("allcause_state_byweek_suppressed_25mar21.xlsx")
icd10categories <- read_excel("icd10categories_updated_321.xlsx")
#weekly state data
all_cause_state_byweek <- allcause_state_byweek %>%
rename(weekly_count_statelevel = count)
all_cause_state_by_week <- all_cause_state_byweek %>%
separate(mmwr_week, into = c("MMWRyear", "MMWRweek"), sep = "-", convert = TRUE) %>%
mutate(
date = MMWRweek2Date(MMWRyear = MMWRyear, MMWRweek = MMWRweek, MMWRday = 1),
statecount_numeric = as.numeric(ifelse(weekly_count_statelevel == "<5", NA, weekly_count_statelevel)),
statecount_suppressed = weekly_count_statelevel == "<5",
category = as.character(category)
)
all_cause_state_by_week %>%
group_by(category) %>%
reframe(missing_total = sum(statecount_suppressed),
count_total = sum(statecount_numeric, na.rm = TRUE))## # A tibble: 10 × 3
## category missing_total count_total
## <chr> <int> <dbl>
## 1 cardiovascular 0 373147
## 2 covid 1 58262
## 3 digestive 0 182701
## 4 infectious 0 112279
## 5 infectious respiratory 0 29816
## 6 injury 0 331205
## 7 other 0 1740863
## 8 pneumonia 0 97954
## 9 renal 0 102199
## 10 respiratory 0 164424all_cause_state_by_week## # A tibble: 3,140 × 7
## MMWRyear MMWRweek category weekly_count_statelevel date
## <int> <int> <chr> <chr> <date>
## 1 2019 1 cardiovascular 843 2018-12-30
## 2 2019 1 digestive 441 2018-12-30
## 3 2019 1 infectious 299 2018-12-30
## 4 2019 1 infectious respiratory 240 2018-12-30
## 5 2019 1 injury 691 2018-12-30
## 6 2019 1 other 4090 2018-12-30
## 7 2019 1 pneumonia 427 2018-12-30
## 8 2019 1 renal 220 2018-12-30
## 9 2019 1 respiratory 522 2018-12-30
## 10 2019 1 covid 0 2018-12-30
## # ℹ 3,130 more rows
## # ℹ 2 more variables: statecount_numeric <dbl>, statecount_suppressed <lgl>str(icd10categories)## tibble [869 × 5] (S3: tbl_df/tbl/data.frame)
## $ icd10code : chr [1:869] "A00" "A01" "A02" "A03" ...
## $ category_original : chr [1:869] "infectious" "infectious" "infectious" "infectious" ...
## $ category_1113update: chr [1:869] "infectious" "infectious" "infectious" "infectious" ...
## $ category_123update : chr [1:869] "infectious" "infectious" "infectious" "infectious" ...
## $ category_0321update: chr [1:869] "infectious" "infectious" "infectious" "infectious" ...names(icd10categories)## [1] "icd10code" "category_original" "category_1113update"
## [4] "category_123update" "category_0321update"dim(icd10categories)## [1] 869 5summary(icd10categories)## icd10code category_original category_1113update category_123update
## Length:869 Length:869 Length:869 Length:869
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
## category_0321update
## Length:869
## Class :character
## Mode :charactersapply(icd10categories, class)## icd10code category_original category_1113update category_123update
## "character" "character" "character" "character"
## category_0321update
## "character"glimpse(icd10categories)## Rows: 869
## Columns: 5
## $ icd10code <chr> "A00", "A01", "A02", "A03", "A04", "A05", "A06", "…
## $ category_original <chr> "infectious", "infectious", "infectious", "infecti…
## $ category_1113update <chr> "infectious", "infectious", "infectious", "infecti…
## $ category_123update <chr> "infectious", "infectious", "infectious", "infecti…
## $ category_0321update <chr> "infectious", "infectious", "infectious", "infecti…colSums(is.na(icd10categories))## icd10code category_original category_1113update category_123update
## 0 418 418 1
## category_0321update
## 0str(allcause_state_byweek)## tibble [3,140 × 3] (S3: tbl_df/tbl/data.frame)
## $ mmwr_week: chr [1:3140] "2019-01" "2019-01" "2019-01" "2019-01" ...
## $ category : chr [1:3140] "cardiovascular" "digestive" "infectious" "infectious respiratory" ...
## $ count : chr [1:3140] "843" "441" "299" "240" ...names(allcause_state_byweek)## [1] "mmwr_week" "category" "count"dim(allcause_state_byweek)## [1] 3140 3summary(allcause_state_byweek)## mmwr_week category count
## Length:3140 Length:3140 Length:3140
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :charactersapply(allcause_state_byweek,class)## mmwr_week category count
## "character" "character" "character"glimpse(allcause_state_byweek)## Rows: 3,140
## Columns: 3
## $ mmwr_week <chr> "2019-01", "2019-01", "2019-01", "2019-01", "2019-01", "2019…
## $ category <chr> "cardiovascular", "digestive", "infectious", "infectious res…
## $ count <chr> "843", "441", "299", "240", "691", "4090", "427", "220", "52…colSums(is.na(allcause_state_byweek))## mmwr_week category count
## 0 0 0Deficit Function Used
deficit_cumulative <- function(fit, start, end){
if (!"curve_fit" %in% attr(fit, "type"))
stop("This is not the correct deficit_model fit, needs curve fit.")
ind <- which(fit$date %in% seq(start, end, by = "day"))
n <- length(ind) # returns the # of elements in the ind object
A <- matrix(1, n, n)
A[upper.tri(A)] <- 0 # It sets all the upper triangular elements of matrix A (excluding the diagonal) to 0
A <- sweep(A, 2, fit$expected[ind], FUN = "*") # This line multiplies each column (2 means columns) of the matrix A by a corresponding value from fit$expected[ind]
# So, now it's scaling each column of A by expected values from the model fit (fit$expected[ind])
fhat <- matrix(fit$fitted[ind], ncol = 1) # creates a column vector (matrix with 1 column) called fhat from the values in fit$fitted[ind]
# Below, this code is constructing a data frame with:
fit_deficit <- A %*% fhat # multiply design matrix A by fitted values fhat to get cumulative modeled effect.
obs_deficit <- cumsum(fit$observed[ind] - fit$expected[ind])
varcov <- fit$x[ind,] %*% fit$betacov %*% t(fit$x[ind,]) # Compute variance-covariance matrix of the linear predictor based on model design matrix x and estimated betacov
diag(varcov) <- fit$se[ind]^2 # Replace diagonal entries of varcov with squared standard errors — probably to correct for over/under-dispersion or model misspecification.
fitted_se <- sqrt(diag(A %*% varcov %*% t(A)))
sd <- sqrt(diag(A %*% (fit$cov[ind,ind] + diag(fit$log_expected_se[ind]^2)) %*% t(A))) # this computes the variance of the observed cumulative deficit.
# It adds fit$cov (covariance of observed data) and extra uncertainty from fit$log_expected_se.
data.frame(date = fit$date[ind], # The time points
observed = obs_deficit, # Cumulative observed - expected difference
sd = sd, # Modeled cumulative deficit/excess from the model.
fitted = fit_deficit, # Standard deviation of the observed cumulative deficit.
se = fitted_se) # Standard error of the fitted cumulative deficit.
}Deficit Model fits for the different categories
Cardio
exclude_dates <- seq(lubridate::make_date(2020, 3, 1),
lubridate::make_date(2021, 2, 20),
by = "day")
# cardiovascular data
cardio <- all_cause_state_by_week %>%
filter(category == "cardiovascular") %>%
filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
outcome = statecount_numeric)
# deficit model with COVID-19 dates excluded
cardio_deficit <- compute_expected(cardio,
exclude = exclude_dates,
frequency = 52)## Overall death rate is 8.85.cardio_deficit## # A tibble: 312 × 12
## MMWRyear MMWRweek category weekly_count_statelevel date
## <int> <int> <chr> <chr> <date>
## 1 2019 2 cardiovascular 1233 2019-01-06
## 2 2019 3 cardiovascular 1160 2019-01-13
## 3 2019 4 cardiovascular 1201 2019-01-20
## 4 2019 5 cardiovascular 1157 2019-01-27
## 5 2019 6 cardiovascular 1286 2019-02-03
## 6 2019 7 cardiovascular 1261 2019-02-10
## 7 2019 8 cardiovascular 1141 2019-02-17
## 8 2019 9 cardiovascular 1156 2019-02-24
## 9 2019 10 cardiovascular 1208 2019-03-03
## 10 2019 11 cardiovascular 1232 2019-03-10
## # ℹ 302 more rows
## # ℹ 7 more variables: statecount_numeric <dbl>, statecount_suppressed <lgl>,
## # population <dbl>, outcome <dbl>, log_expected_se <dbl>, expected <dbl>,
## # excluded <lgl>cardio_deficit <- cardio_deficit %>%
mutate(date = as.Date(paste(MMWRyear, MMWRweek, 1, sep = "-"), "%Y-%U-%u"))## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `date = as.Date(paste(MMWRyear, MMWRweek, 1, sep = "-"),
## "%Y-%U-%u")`.
## Caused by warning in `strptime()`:
## ! (0-based) yday 369 in year 2020 is invalidcombined_plot_1 <- ggplot(cardio_deficit, aes(x = date)) +
geom_line(aes(y = expected), color = "blue", size = 1.5) + # Expected outcome line
geom_point(aes(y = outcome, color = factor(excluded)), size = 2) + # Observed outcomes
scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
labels = c("Included in Model", "Excluded from Model")) +
labs(title = "Weekly Cardiovascular Hospitalizations in MA",
x = "Date",
y = "Outcome",
color = "Data Point Status") +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.combined_plot_1## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
fit_x <- cardio_deficit %>%
deficit_model(exclude = exclude_dates,
start = as.Date("2019-01-14"),
end = as.Date("2024-12-31"),
knots.per.year = 12,
verbose = FALSE)
# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_x <- excess_plot(fit_x,
title = "Deviations from Expected Counts for Cardio in MA"
)
deficit_plot_x
fit_x$detected_intervals## start end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-09-23 2019-12-02 8.144139 8.669496 0.07650853 12062
## 2 2020-01-20 2020-02-03 8.276251 8.703193 0.14678711 3343
## 3 2020-03-02 2021-02-15 7.265876 8.874462 0.03594969 49893
## 4 2021-12-27 2022-01-31 8.304721 9.024090 0.10569035 6709
## 5 2023-05-01 2023-05-01 9.677494 9.573722 0.26665497 1303
## expected deficit sd
## 1 12840.088 -778.0880 113.31411
## 2 3515.453 -172.4533 59.29126
## 3 60938.774 -11045.7736 246.85780
## 4 7290.145 -581.1446 85.38234
## 5 1289.028 13.9721 35.90303# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit <- deficit_cumulative(fit_x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_1 <- cumulative_deficit %>%
ggplot(aes(date)) +
geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
geom_line(aes(y = observed),
color = "white",
size = 1) +
geom_line(aes(y = observed)) +
geom_point(aes(y = observed)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Date",
y = "Cumulative Deficit Hospitilization-Cardio",
title = "Cumulative Deficit Hospitilization for Cardio Cases in MA",
subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_1
library(patchwork)
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")
combined_plot_1_adj <- combined_plot_1 +
theme(legend.position = "bottom",
legend.justification = c(1, 0),
legend.box.just = "right") +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
deficit_plot_x_adj <- deficit_plot_x +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.cumulative_deficit_plot_1_adj <- cumulative_deficit_plot_1 +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
final_combined_plot_1 <- combined_plot_1_adj /
deficit_plot_x_adj /
cumulative_deficit_plot_1_adj +
plot_layout(ncol = 1, heights = c(1, 1, 1))
final_combined_plot_1## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("Cardio_combined_figure.png",
plot = final_combined_plot_1,
width = 11,
height = 10.5,
dpi = 300,
bg = "white")## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).ggsave("cardio_combined_figure.png",
plot = final_combined_plot_1,
width = 11, # good for standard document width
height = 10.5, # leaves space for three stacked plots
dpi = 300,
bg = "white")## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).Digestive
# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1),
lubridate::make_date(2021, 2, 27),
by = "day")
# digestive data
digestive <- all_cause_state_by_week %>%
filter(category == "digestive") %>%
filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
outcome = statecount_numeric)
# excess model with COVID-19 dates excluded
digestive_deficit <- compute_expected(digestive,
exclude = exclude_dates,
frequency = 52)## Overall death rate is 4.33.digestive_deficit## # A tibble: 312 × 12
## MMWRyear MMWRweek category weekly_count_statelevel date
## <int> <int> <chr> <chr> <date>
## 1 2019 2 digestive 535 2019-01-06
## 2 2019 3 digestive 570 2019-01-13
## 3 2019 4 digestive 553 2019-01-20
## 4 2019 5 digestive 576 2019-01-27
## 5 2019 6 digestive 609 2019-02-03
## 6 2019 7 digestive 599 2019-02-10
## 7 2019 8 digestive 602 2019-02-17
## 8 2019 9 digestive 562 2019-02-24
## 9 2019 10 digestive 550 2019-03-03
## 10 2019 11 digestive 622 2019-03-10
## # ℹ 302 more rows
## # ℹ 7 more variables: statecount_numeric <dbl>, statecount_suppressed <lgl>,
## # population <dbl>, outcome <dbl>, log_expected_se <dbl>, expected <dbl>,
## # excluded <lgl># Plotting
combined_plot_2 <- ggplot(digestive_deficit, aes(x = date)) +
geom_line(aes(y = expected), color = "blue", size = 1.5) + # Expected outcome line
geom_point(aes(y = outcome, color = factor(excluded)), size = 2) + # Observed outcomes
scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
labels = c("Included in Model", "Excluded from Model")) +
labs(title = "Weekly digestive ED Hospitilizations in MA",
x = "Date",
y = "Outcome",
color = "Data Point Status") +
theme_minimal()
combined_plot_2
fit_2x <- digestive_deficit %>%
deficit_model(exclude = exclude_dates,
start = min(.$date),
end = max(.$date),
knots.per.year = 12,
verbose = FALSE)
# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_2x <- excess_plot(fit_2x,
title = "Deviations from Expected Counts for Digestive in MA"
)
deficit_plot_2x
fit_2x$detected_intervals## start end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-09-15 2019-10-13 3.986860 4.261767 0.07956446 2684
## 2 2020-03-01 2020-09-06 3.519113 4.442692 0.03432838 13267
## 3 2020-10-11 2021-03-14 3.714835 4.193817 0.03680019 11504
## 4 2021-12-05 2022-01-30 3.759757 4.186787 0.05877984 4556
## 5 2022-05-01 2022-05-22 4.368984 4.646495 0.09288423 2353
## 6 2022-09-18 2022-10-02 4.263148 4.543836 0.10606203 1722
## expected deficit sd
## 1 2869.071 -185.0706 53.56371
## 2 16748.877 -3481.8769 129.41745
## 3 12987.299 -1483.2986 113.96183
## 4 5073.467 -517.4671 71.22827
## 5 2502.459 -149.4588 50.02458
## 6 1835.377 -113.3774 42.84131# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit <- deficit_cumulative(fit_2x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_2 <- cumulative_deficit %>%
ggplot(aes(date)) +
geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
geom_line(aes(y = observed),
color = "white",
size = 1) +
geom_line(aes(y = observed)) +
geom_point(aes(y = observed)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Date",
y = "Cumulative Deficit Hospitilization-Digestive",
title = "Cumulative Deficit Hospitilization for Digestive Cases in MA",
subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_2
library(patchwork)
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")
combined_plot_2_adj <- combined_plot_2 +
theme(legend.position = "bottom",
legend.justification = c(1, 0),
legend.box.just = "right") +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
deficit_plot_2x_adj <- deficit_plot_2x +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.cumulative_deficit_plot_2_adj <- cumulative_deficit_plot_2 +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
final_combined_plot_2 <- combined_plot_2_adj /
deficit_plot_2x_adj /
cumulative_deficit_plot_2_adj +
plot_layout(ncol = 1, heights = c(1, 1, 1))
final_combined_plot_2
ggsave("Digestive_combined_figure.png",
plot = final_combined_plot_2,
width = 11,
height = 10.5,
dpi = 300,
bg = "white")Infectious
# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1),
lubridate::make_date(2021, 6, 5),
by = "day")
# infectious data
infectious <- all_cause_state_by_week %>%
filter(category == "infectious") %>%
filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
outcome = statecount_numeric)
# deficit model with COVID-19 dates excluded
infectious_deficit <- compute_expected(infectious,
exclude = exclude_dates,
frequency = 52)## Overall death rate is 2.66.# Plotting
combined_plot_3 <- ggplot(infectious_deficit, aes(x = date)) +
geom_line(aes(y = expected), color = "blue", size = 1.5) + # Expected outcome line
geom_point(aes(y = outcome, color = factor(excluded)), size = 2) + # Observed outcomes
scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
labels = c("Included in Model", "Excluded from Model")) +
labs(title = "Weekly Infectious ED Hospitilizations in MA",
x = "Date",
y = "Outcome",
color = "Data Point Status") +
theme_minimal()
combined_plot_3
fit_3x <- infectious_deficit %>%
deficit_model(exclude = exclude_dates,
start = as.Date("2020-03-01"),
end = as.Date("2024-12-31"),
knots.per.year = 12,
verbose = FALSE)
deficit_plot_3x <- excess_plot(fit_3x,
title = "Deviations from Expected Counts for Infectious in MA"
)
deficit_plot_3x
fit_3x$detected_intervals## start end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2020-05-03 2020-08-02 2.384625 2.718529 0.03797626 4495
## 2 2020-10-04 2021-05-30 2.231734 2.726452 0.02405327 10517
## 3 2022-01-23 2022-02-27 2.530161 2.785967 0.05872481 2044
## 4 2022-12-25 2023-01-22 2.565316 2.820316 0.06472516 1727
## 5 2023-06-18 2023-07-02 2.542539 2.771474 0.08283312 1027
## 6 2024-10-20 2024-12-15 2.525210 2.768028 0.04779399 3060
## expected deficit sd
## 1 5124.405 -629.40467 71.58495
## 2 12848.349 -2331.34905 113.35056
## 3 2250.654 -206.65422 47.44106
## 4 1898.669 -171.66889 43.57372
## 5 1119.473 -92.47266 33.45852
## 6 3354.243 -294.24302 57.91583# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit <- deficit_cumulative(fit_3x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_3 <- cumulative_deficit %>%
ggplot(aes(date)) +
geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
geom_line(aes(y = observed),
color = "white",
size = 1) +
geom_line(aes(y = observed)) +
geom_point(aes(y = observed)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Date",
y = "Cumulative Deficit Hospitilization-Infectious",
title = "Cumulative Deficit Hospitilization for Infectious Cases in MA",
subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_3
library(patchwork)
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")
combined_plot_3_adj <- combined_plot_3 +
theme(legend.position = "bottom",
legend.justification = c(1, 0),
legend.box.just = "right") +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
deficit_plot_3x_adj <- deficit_plot_3x +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.cumulative_deficit_plot_3_adj <- cumulative_deficit_plot_3 +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
final_combined_plot_3 <- combined_plot_3_adj /
deficit_plot_3x_adj /
cumulative_deficit_plot_3_adj +
plot_layout(ncol = 1, heights = c(1, 1, 1))
final_combined_plot_3
ggsave("Infectious_combined_figure.png",
plot = final_combined_plot_2,
width = 11,
height = 10.5,
dpi = 300,
bg = "white")infectious_respiratory
# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1),
lubridate::make_date(2021, 6, 12),
by = "day")
# infectious respiratory data
infectious_respiratory <- all_cause_state_by_week %>%
filter(category == "infectious respiratory") %>%
filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
outcome = statecount_numeric)
# excess model with COVID-19 dates excluded
infectious_respiratory_deficit <- compute_expected(infectious_respiratory,
exclude = exclude_dates,
frequency = 52)## Overall death rate is 0.699.# Plotting
combined_plot_4 <- ggplot(infectious_respiratory_deficit, aes(x = date)) +
geom_line(aes(y = expected), color = "blue", size = 1.5) + # Expected outcome line
geom_point(aes(y = outcome, color = factor(excluded)), size = 2) + # Observed outcomes
scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
labels = c("Included in Model", "Excluded from Model")) +
labs(title = "Weekly Infectious Respiratory ED Hospitilizations in MA",
x = "Date",
y = "Outcome",
color = "Data Point Status") +
theme_minimal()
combined_plot_4
# -- Fitting excess model
fit_4x <- infectious_respiratory_deficit %>%
deficit_model(exclude = exclude_dates,
start = as.Date("2020-03-01"),
end = as.Date("2025-12-31"),
knots.per.year = 12,
verbose = FALSE)
# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_4x <- excess_plot(fit_4x,
title = "Deviations from Expected Counts for Infectious Respiratory in MA"
)
deficit_plot_4x
# -- Intervals of inordinate mortality found by the excess model
fit_4x$detected_intervals## start end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2020-03-29 2021-05-30 0.1241042 0.8568669 0.01013141 1036
## 2 2021-10-31 2022-03-20 0.4714432 1.3271874 0.02166535 1333
## 3 2023-01-15 2023-02-26 0.6535837 1.2568241 0.03651720 616
## 4 2023-04-09 2023-04-23 0.4208683 0.7014588 0.04167251 170
## 5 2024-11-03 2024-11-10 0.6424430 1.0037101 0.06105179 173
## expected deficit sd
## 1 7152.9723 -6116.97228 84.57525
## 2 3752.6065 -2419.60646 61.25852
## 3 1184.5517 -568.55175 34.41732
## 4 283.3380 -113.33805 16.83265
## 5 270.2837 -97.28366 16.44031# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit <- deficit_cumulative(fit_4x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_4 <- cumulative_deficit %>%
ggplot(aes(date)) +
geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
geom_line(aes(y = observed),
color = "white",
size = 1) +
geom_line(aes(y = observed)) +
geom_point(aes(y = observed)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Date",
y = "Cumulative Deficit Hospitilization-infectious respiratory",
title = "Cumulative Deficit Hospitilization for infectious respiratory Cases in MA",
subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_4
library(patchwork)
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")
combined_plot_4_adj <- combined_plot_4 +
theme(legend.position = "bottom",
legend.justification = c(1, 0),
legend.box.just = "right") +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
deficit_plot_4x_adj <- deficit_plot_4x +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.cumulative_deficit_plot_4_adj <- cumulative_deficit_plot_4 +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
final_combined_plot_4 <- combined_plot_4_adj /
deficit_plot_4x_adj /
cumulative_deficit_plot_4_adj +
plot_layout(ncol = 1, heights = c(1, 1, 1))
final_combined_plot_4
ggsave("Infectious Respiratory_combined_figure.png",
plot = final_combined_plot_4,
width = 11, # good for standard document width
height = 10.5, # leaves space for three stacked plots
dpi = 300,
bg = "white")injury
# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1),
lubridate::make_date(2021, 4, 10),
by = "day")
# injury data
injury <- all_cause_state_by_week %>%
filter(category == "injury") %>%
filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
outcome = statecount_numeric)
# excess model with COVID-19 dates excluded
injury_deficit <- compute_expected(injury,
exclude = exclude_dates,
frequency = 52)## Overall death rate is 7.86.# Plotting
combined_plot_5 <- ggplot(injury_deficit, aes(x = date)) +
geom_line(aes(y = expected), color = "blue", size = 1.5) + # Expected outcome line
geom_point(aes(y = outcome, color = factor(excluded)), size = 2) + # Observed outcomes
scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
labels = c("Included in Model", "Excluded from Model")) +
labs(title = "Weekly Injury ED Hospitilizations in MA",
x = "Date",
y = "Outcome",
color = "Data Point Status") +
theme_minimal()
combined_plot_5
# -- Fitting excess model
fit_5x <- injury_deficit %>%
deficit_model(exclude = exclude_dates,
start = as.Date("2020-03-01"),
end = as.Date("2024-12-31"),
knots.per.year = 12,
verbose = FALSE)
# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_5x <- excess_plot(fit_5x,
title = "Deviations from Expected Counts for Injuries in MA"
)
deficit_plot_5x
# -- Intervals of inordinate mortality found by the excess model
fit_5x$detected_intervals## start end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2020-03-08 2020-09-20 6.347086 7.766319 0.04459825 24783
## 2 2020-10-04 2021-02-14 6.613078 7.488809 0.05273521 17808
## 3 2021-03-07 2021-04-11 7.120101 7.580370 0.09686766 5752
## 4 2021-05-09 2021-05-16 7.731598 7.976285 0.17210544 2082
## 5 2021-07-25 2021-08-08 7.979167 8.372155 0.14396843 3223
## 6 2021-12-19 2022-01-16 7.140402 7.682093 0.10682282 4807
## 7 2022-07-31 2022-08-21 8.268205 8.668967 0.12687117 4453
## expected deficit sd
## 1 30324.572 -5541.57214 174.13952
## 2 20166.209 -2358.20891 142.00778
## 3 6123.830 -371.83007 78.25490
## 4 2147.891 -65.89052 46.34534
## 5 3381.738 -158.73836 58.15272
## 6 5171.673 -364.67274 71.91434
## 7 4668.838 -215.83814 68.32890# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit <- deficit_cumulative(fit_5x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_5 <- cumulative_deficit %>%
ggplot(aes(date)) +
geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
geom_line(aes(y = observed),
color = "white",
size = 1) +
geom_line(aes(y = observed)) +
geom_point(aes(y = observed)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Date",
y = "Cumulative Deficit Hospitilization-Injury",
title = "Cumulative Deficit Hospitilization for Injury Cases in MA",
subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_5
library(patchwork)
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")
combined_plot_5_adj <- combined_plot_5 +
theme(legend.position = "bottom",
legend.justification = c(1, 0),
legend.box.just = "right") +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
deficit_plot_5x_adj <- deficit_plot_5x +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.cumulative_deficit_plot_5_adj <- cumulative_deficit_plot_5 +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
final_combined_plot_5 <- combined_plot_5_adj /
deficit_plot_5x_adj /
cumulative_deficit_plot_5_adj +
plot_layout(ncol = 1, heights = c(1, 1, 1))
final_combined_plot_5
ggsave("Injury_combined_figure.png",
plot = final_combined_plot_5,
width = 11, # good for standard document width
height = 10.5, # leaves space for three stacked plots
dpi = 300,
bg = "white")pneumonia
# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1),
lubridate::make_date(2022, 5, 7),
by = "day")
# pneumonia data
pneumonia <- all_cause_state_by_week %>%
filter(category == "pneumonia") %>%
filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
outcome = statecount_numeric)
# excess model with COVID-19 dates excluded
pneumonia_deficit <- compute_expected(pneumonia,
exclude = exclude_dates,
frequency = 52)## Overall death rate is 2.32.# Plotting
combined_plot_6 <- ggplot(pneumonia_deficit, aes(x = date)) +
geom_line(aes(y = expected), color = "blue", size = 1.5) + # Expected outcome line
geom_point(aes(y = outcome, color = factor(excluded)), size = 2) + # Observed outcomes
scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
labels = c("Included in Model", "Excluded from Model")) +
labs(title = "Weekly Pneumonia ED Hospitilizations in MA",
x = "Date",
y = "Outcome",
color = "Data Point Status") +
theme_minimal()
combined_plot_6
# -- Fitting excess model
fit_6x <- pneumonia_deficit %>%
deficit_model(exclude = exclude_dates,
start = as.Date("2020-03-01"),
end = as.Date("2024-12-31"),
knots.per.year = 12,
verbose = FALSE)
# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_6x <- excess_plot(fit_6x,
title = "Deviations from Expected Counts for Pneumonia in MA"
)
deficit_plot_6x
# -- Intervals of inordinate mortality found by the excess model
fit_6x$detected_intervals## start end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2020-04-26 2021-07-04 1.489661 2.618500 0.01756974 12636
## 2 2021-09-05 2022-05-22 1.947460 2.813849 0.02345135 9964
## 3 2022-07-17 2022-08-28 1.694437 2.049678 0.04663405 1597
## 4 2022-10-30 2022-11-13 2.562345 2.880032 0.08443981 1035
## 5 2023-01-22 2023-02-19 2.407862 2.978610 0.06651675 1621
## 6 2023-04-23 2023-05-21 2.294970 2.662542 0.06288869 1545
## 7 2024-01-28 2024-02-11 2.619286 2.974941 0.08581986 1058
## expected deficit sd
## 1 22211.329 -9575.3287 149.03466
## 2 14396.797 -4432.7975 119.98666
## 3 1931.814 -334.8138 43.95240
## 4 1163.322 -128.3221 34.10751
## 5 2005.234 -384.2341 44.77984
## 6 1792.454 -247.4540 42.33738
## 7 1201.659 -143.6586 34.66495# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit <- deficit_cumulative(fit_6x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_6 <- cumulative_deficit %>%
ggplot(aes(date)) +
geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
geom_line(aes(y = observed),
color = "white",
size = 1) +
geom_line(aes(y = observed)) +
geom_point(aes(y = observed)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Date",
y = "Cumulative Deficit Hospitilization-pneumonia",
title = "Cumulative Deficit Hospitilization for pneumonia Cases in MA",
subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_6
library(patchwork)
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")
combined_plot_6_adj <- combined_plot_6 +
theme(legend.position = "bottom",
legend.justification = c(1, 0),
legend.box.just = "right") +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
deficit_plot_6x_adj <- deficit_plot_6x +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.cumulative_deficit_plot_6_adj <- cumulative_deficit_plot_6 +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
final_combined_plot_6 <- combined_plot_6_adj /
deficit_plot_6x_adj /
cumulative_deficit_plot_6_adj +
plot_layout(ncol = 1, heights = c(1, 1, 1))
final_combined_plot_6
ggsave("Pneumonia_combined_figure.png",
plot = final_combined_plot_6,
width = 11,
height = 10.5,
dpi = 300,
bg = "white")Renal
# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1),
lubridate::make_date(2021, 2, 27),
by = "day")
# renal data
renal <- all_cause_state_by_week %>%
filter(category == "renal") %>%
filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
outcome = statecount_numeric)
# excess model with COVID-19 dates excluded
renal_deficit <- compute_expected(renal,
exclude = exclude_dates,
frequency = 52)## Overall death rate is 2.42.# Plotting
combined_plot_7 <- ggplot(renal_deficit, aes(x = date)) +
geom_line(aes(y = expected), color = "blue", size = 1.5) + # Expected outcome line
geom_point(aes(y = outcome, color = factor(excluded)), size = 2) + # Observed outcomes
scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
labels = c("Included in Model", "Excluded from Model")) +
labs(title = "Weekly Renal ED Hospitilizations in MA",
x = "Date",
y = "Outcome",
color = "Data Point Status") +
theme_minimal()
combined_plot_7
# -- Fitting excess model
fit_7x <- renal_deficit %>%
deficit_model(exclude = exclude_dates,
start = as.Date("2020-03-01"),
end = as.Date("2024-12-31"),
knots.per.year = 12,
verbose = FALSE)
# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_7x <- excess_plot(fit_7x,
title = "Deviations from Expected Counts for Renal in MA"
)
deficit_plot_7x
# -- Intervals of inordinate mortality found by the excess model
fit_7x$detected_intervals## start end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2020-03-08 2020-07-12 1.778983 2.317654 0.03009932 4551
## 2 2020-09-06 2020-09-27 2.194704 2.341201 0.06593235 1182
## 3 2020-11-01 2020-12-20 1.925472 2.186413 0.04505369 2074
## 4 2021-01-03 2021-02-28 1.965702 2.247830 0.04306949 2382
## 5 2022-01-02 2022-01-02 2.124147 2.335657 0.13170849 286
## 6 2022-05-15 2022-05-29 2.364289 2.607537 0.08034593 955
## expected deficit sd
## 1 5929.0293 -1378.02930 77.00019
## 2 1260.8985 -78.89848 35.50913
## 3 2355.0697 -281.06966 48.52906
## 4 2723.8769 -341.87691 52.19077
## 5 314.4782 -28.47817 17.73353
## 6 1053.2544 -98.25443 32.45388# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit <- deficit_cumulative(fit_7x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_7 <- cumulative_deficit %>%
ggplot(aes(date)) +
geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
geom_line(aes(y = observed),
color = "white",
size = 1) +
geom_line(aes(y = observed)) +
geom_point(aes(y = observed)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Date",
y = "Cumulative Deficit Hospitilization-renal",
title = "Cumulative Deficit Hospitilization for renal Cases in MA",
subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_7
library(patchwork)
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")
combined_plot_7_adj <- combined_plot_7 +
theme(legend.position = "bottom",
legend.justification = c(1, 0),
legend.box.just = "right") +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
deficit_plot_7x_adj <- deficit_plot_7x +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.cumulative_deficit_plot_7_adj <- cumulative_deficit_plot_7 +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
final_combined_plot_7 <- combined_plot_7_adj /
deficit_plot_7x_adj /
cumulative_deficit_plot_7_adj +
plot_layout(ncol = 1, heights = c(1, 1, 1))
final_combined_plot_7
ggsave("Renal_combined_figure.png",
plot = final_combined_plot_7,
width = 11, # good for standard document width
height = 10.5, # leaves space for three stacked plots
dpi = 300,
bg = "white")Respiratory
# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1),
lubridate::make_date(2021, 6, 12),
by = "day")
# respiratory data
respiratory <- all_cause_state_by_week %>%
filter(category == "respiratory") %>%
filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
outcome = statecount_numeric)
# excess model with COVID-19 dates excluded
respiratory_deficit <- compute_expected(respiratory,
exclude = exclude_dates,
frequency = 52)## Overall death rate is 3.89.# Plotting
combined_plot_8 <- ggplot(respiratory_deficit, aes(x = date)) +
geom_line(aes(y = expected), color = "blue", size = 1.5) + # Expected outcome line
geom_point(aes(y = outcome, color = factor(excluded)), size = 2) + # Observed outcomes
scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
labels = c("Included in Model", "Excluded from Model")) +
labs(title = "Weekly Respiratory ED Hospitilizations in MA",
x = "Date",
y = "Outcome",
color = "Data Point Status") +
theme_minimal()
combined_plot_8
# -- Fitting excess model
fit_8x <- respiratory_deficit %>%
deficit_model(exclude = exclude_dates,
start = as.Date("2020-03-01"),
end = as.Date("2024-12-31"),
knots.per.year = 12,
verbose = FALSE)
# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_8x <- excess_plot(fit_8x,
title = "Deviations from Expected Counts for Respiratory in MA"
)
deficit_plot_8x
# -- Intervals of inordinate mortality found by the excess model
fit_8x$detected_intervals## start end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2020-03-08 2021-05-30 2.675694 4.168126 0.02182343 23417
## 2 2021-11-21 2022-03-27 3.632627 4.441859 0.04166917 9293
## 3 2023-04-23 2023-05-07 4.186401 4.696742 0.10783183 1691
## expected deficit sd
## 1 36478.39 -13061.3943 190.99318
## 2 11363.18 -2070.1786 106.59821
## 3 1897.14 -206.1403 43.55617# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit <- deficit_cumulative(fit_8x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_8 <- cumulative_deficit %>%
ggplot(aes(date)) +
geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
geom_line(aes(y = observed),
color = "white",
size = 1) +
geom_line(aes(y = observed)) +
geom_point(aes(y = observed)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Date",
y = "Cumulative Deficit Hospitilization-respiratory",
title = "Cumulative Deficit Hospitilization for respiratory Cases in MA",
subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_8
library(patchwork)
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")
combined_plot_8_adj <- combined_plot_8 +
theme(legend.position = "bottom",
legend.justification = c(1, 0),
legend.box.just = "right") +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
deficit_plot_8x_adj <- deficit_plot_8x +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.cumulative_deficit_plot_8_adj <- cumulative_deficit_plot_8 +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
final_combined_plot_8 <- combined_plot_8_adj /
deficit_plot_8x_adj /
cumulative_deficit_plot_8_adj +
plot_layout(ncol = 1, heights = c(1, 1, 1))
final_combined_plot_8
ggsave("Respiratory_combined_figure.png",
plot = final_combined_plot_8,
width = 11,
height = 10.5,
dpi = 300,
bg = "white")other
# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1),
lubridate::make_date(2021, 6, 5),
by = "day")
# other data
other <- all_cause_state_by_week %>%
filter(category == "other") %>%
filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
outcome = statecount_numeric)
# excess model with COVID-19 dates excluded
other_deficit <- compute_expected(other,
exclude = exclude_dates,
frequency = 52)## Overall death rate is 41.3.# Plotting
combined_plot_9 <- ggplot(other_deficit, aes(x = date)) +
geom_line(aes(y = expected), color = "blue", size = 1.5) + # Expected outcome line
geom_point(aes(y = outcome, color = factor(excluded)), size = 2) + # Observed outcomes
scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
labels = c("Included in Model", "Excluded from Model")) +
labs(title = "Weekly Other ED Hospitilizations in MA",
x = "Date",
y = "Outcome",
color = "Data Point Status") +
theme_minimal()
combined_plot_9
# -- Fitting excess model
fit_9x <- other_deficit %>%
deficit_model(exclude = exclude_dates,
start = as.Date("2020-03-01"),
end = as.Date("2024-12-31"),
knots.per.year = 12,
verbose = FALSE)
# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_9x <- excess_plot(fit_9x,
title = "Deviations from Expected Counts for Other in MA"
)
deficit_plot_9x
# -- Intervals of inordinate mortality found by the excess model
fit_9x$detected_intervals## start end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2020-03-08 2021-05-23 35.63609 41.62600 0.06950269 307080
## 2 2021-08-22 2021-10-10 40.63174 42.79761 0.19933061 43766
## 3 2021-11-21 2022-02-13 36.00309 40.26895 0.15167805 63018
## 4 2022-03-06 2022-04-10 41.29956 43.11281 0.23101320 33364
## 5 2022-04-24 2022-05-08 42.59682 44.33594 0.33130392 17206
## expected deficit sd
## 1 358695.65 -51615.6461 598.9121
## 2 46098.94 -2332.9402 214.7066
## 3 70484.75 -7466.7499 265.4896
## 4 34828.84 -1464.8444 186.6249
## 5 17908.48 -702.4761 133.8226# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit <- deficit_cumulative(fit_9x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_9 <- cumulative_deficit %>%
ggplot(aes(date)) +
geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
geom_line(aes(y = observed),
color = "white",
size = 1) +
geom_line(aes(y = observed)) +
geom_point(aes(y = observed)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Date",
y = "Cumulative Deficit Hospitilization-other",
title = "Cumulative Deficit Hospitilization for other Cases in MA",
subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_9
library(patchwork)
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")
combined_plot_9_adj <- combined_plot_9 +
theme(legend.position = "bottom",
legend.justification = c(1, 0),
legend.box.just = "right") +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
deficit_plot_9x_adj <- deficit_plot_9x +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.cumulative_deficit_plot_9_adj <- cumulative_deficit_plot_9 +
labs(x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
final_combined_plot_9 <- combined_plot_9_adj /
deficit_plot_9x_adj /
cumulative_deficit_plot_9_adj +
plot_layout(ncol = 1, heights = c(1, 1, 1))
final_combined_plot_9
ggsave("Other_combined_figure.png",
plot = final_combined_plot_5,
width = 11,
height = 10.5,
dpi = 300,
bg = "white")Different Approuch
library(ggplot2)
library(patchwork)
library(cowplot)
library(scales)
# Common scale settings
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(x_limits[1], x_limits[2], by = "1 year")
category_labels <- c("Cardiovascular", "Digestive", "Infectious", "Infectious Respiratory",
"Injury", "Pneumonia", "Renal", "Respiratory", "Other")
# ----- CLEANERS -----
clean_combined <- function(plots, labels) {
mapply(function(p, label) {
p +
labs(title = label, x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y") +
theme_minimal() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, size = 11),
axis.text.x = element_text(size = 8))
}, plots, labels, SIMPLIFY = FALSE)
}
clean_deficit <- function(plot_list, labels) {
mapply(function(p, label) {
p +
labs(title = label, x = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, size = 11),
axis.text.x = element_text(size = 8))
}, plot_list, labels, SIMPLIFY = FALSE)
}
clean_cumulative <- function(plots, labels) {
mapply(function(p, label) {
p +
labs(title = label, x = NULL, y = NULL, subtitle = NULL) +
scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y") +
scale_y_continuous(labels = comma) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, size = 11),
axis.text = element_text(size = 8))
}, plots, labels, SIMPLIFY = FALSE)
}
# ----- APPLY CLEANING -----
combined_clean <- clean_combined(list(
combined_plot_1, combined_plot_2, combined_plot_3,
combined_plot_4, combined_plot_5, combined_plot_6,
combined_plot_7, combined_plot_8, combined_plot_9
), category_labels)
deficit_clean <- clean_deficit(list(
deficit_plot_x, deficit_plot_2x, deficit_plot_3x,
deficit_plot_4x, deficit_plot_5x, deficit_plot_6x,
deficit_plot_7x, deficit_plot_8x, deficit_plot_9x
), category_labels)## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.cumulative_clean <- clean_cumulative(list(
cumulative_deficit_plot_1, cumulative_deficit_plot_2, cumulative_deficit_plot_3,
cumulative_deficit_plot_4, cumulative_deficit_plot_5, cumulative_deficit_plot_6,
cumulative_deficit_plot_7, cumulative_deficit_plot_8, cumulative_deficit_plot_9
), category_labels)## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.# ----- LEGEND FOR COMBINED -----
legend_shared <- get_legend(
combined_plot_1 +
theme(legend.position = "bottom",
legend.title = element_text(size = 10),
legend.text = element_text(size = 9))
)## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).## Warning in get_plot_component(plot, "guide-box"): Multiple components found;
## returning the first one. To return all, use `return_all = TRUE`.# ----- FINAL COMBINED PLOT (balanced, readable) -----
combined_final <- wrap_plots(combined_clean, ncol = 3) +
plot_layout(guides = "collect") &
theme(legend.position = "bottom")
combined_final <- combined_final +
plot_annotation(
title = "Weekly Observed vs Expected Hospitalizations by Cause",
theme = theme(plot.title = element_text(hjust = 0.5, size = 15))
)
ggsave("combined_all_categories.png",
plot = combined_final,
width = 12, height = 10, dpi = 300, bg = "white")## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).deficit_final <- wrap_plots(deficit_clean, ncol = 3) +
plot_annotation(
title = "Deviations from Expected Weekly Admissions",
theme = theme(plot.title = element_text(hjust = 0.5, size = 14))
)
ggsave("deficit_all_categories.png",
plot = deficit_final,
width = 11, height = 9.5, dpi = 300, bg = "white")
cumulative_grid <- wrap_plots(cumulative_clean, ncol = 3)
y_label <- ggdraw() +
draw_label("Cumulative Deficit in Hospitalizations", angle = 90, size = 12)
cumulative_final <- plot_grid(
y_label, cumulative_grid,
rel_widths = c(0.04, 0.96),
nrow = 1
)
cumulative_final <- cumulative_final +
plot_annotation(
title = "Cumulative Deficit in Hospital Admissions by Cause (2020–2024)",
theme = theme(plot.title = element_text(hjust = 0.5, size = 15))
)
ggsave("cumulative_all_categories.png",
plot = cumulative_final,
width = 12, height = 10, dpi = 300, bg = "white")
combined_final## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
deficit_final
cumulative_final
ggsave("deficit_plot_x.png", plot = deficit_plot_x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_2x.png", plot = deficit_plot_2x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_3x.png", plot = deficit_plot_3x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_4x.png", plot = deficit_plot_4x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_5x.png", plot = deficit_plot_5x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_6x.png", plot = deficit_plot_6x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_7x.png", plot = deficit_plot_7x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_8x.png", plot = deficit_plot_8x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_9x.png", plot = deficit_plot_9x, width = 8, height = 6, dpi = 300)Cumulative deficit interactive plots
cumulative_deficit <- deficit_cumulative(fit_x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
direct_plotly <- plot_ly(data = cumulative_deficit) %>%
add_ribbons(x = ~date,
ymin = ~observed - 2*sd,
ymax = ~observed + 2*sd,
name = "95% CI",
fillcolor = "rgba(180, 180, 180, 0.5)",
line = list(color = "transparent"),
hoverinfo = "none") %>%
add_lines(x = ~date,
y = ~observed,
name = "Observed",
line = list(color = "black", width = 2),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
add_markers(x = ~date,
y = ~observed,
name = "Data Points",
marker = list(color = "black", size = 4),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
layout(title = list(text = "Cumulative Deficit Hospitalization for Cardio Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
xaxis = list(title = "Date"),
yaxis = list(title = "Cumulative Deficit Hospitalization-Cardio",
tickformat = ","),
hoverlabel = list(bgcolor = "white",
font = list(family = "Arial", size = 12)),
hovermode = "closest")
direct_plotly# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit <- deficit_cumulative(fit_2x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_2 <- cumulative_deficit %>%
ggplot(aes(date)) +
geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
geom_line(aes(y = observed),
color = "white",
size = 1) +
geom_line(aes(y = observed)) +
geom_point(aes(y = observed)) +
scale_y_continuous(labels = scales::comma) +
labs(x = "Date",
y = "Cumulative Deficit Hospitilization-Digestive",
title = "Cumulative Deficit Hospitilization for Digestive Cases in MA",
subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_2
cumulative_deficit <- deficit_cumulative(fit_2x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
direct_plotly <- plot_ly(data = cumulative_deficit) %>%
add_ribbons(x = ~date,
ymin = ~observed - 2*sd,
ymax = ~observed + 2*sd,
name = "95% CI",
fillcolor = "rgba(180, 180, 180, 0.5)",
line = list(color = "transparent"),
hoverinfo = "none") %>%
add_lines(x = ~date,
y = ~observed,
name = "Observed",
line = list(color = "black", width = 2),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
add_markers(x = ~date,
y = ~observed,
name = "Data Points",
marker = list(color = "black", size = 4),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
layout(title = list(text = "Cumulative Deficit Hospitalization for Digestive Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
xaxis = list(title = "Date"),
yaxis = list(title = "Cumulative Deficit Hospitalization-Digestive",
tickformat = ","),
hoverlabel = list(bgcolor = "white",
font = list(family = "Arial", size = 12)),
hovermode = "closest")
direct_plotly#3
cumulative_deficit <- deficit_cumulative(fit_3x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
direct_plotly <- plot_ly(data = cumulative_deficit) %>%
add_ribbons(x = ~date,
ymin = ~observed - 2*sd,
ymax = ~observed + 2*sd,
name = "95% CI",
fillcolor = "rgba(180, 180, 180, 0.5)",
line = list(color = "transparent"),
hoverinfo = "none") %>%
add_lines(x = ~date,
y = ~observed,
name = "Observed",
line = list(color = "black", width = 2),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
add_markers(x = ~date,
y = ~observed,
name = "Data Points",
marker = list(color = "black", size = 4),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
layout(title = list(text = "Cumulative Deficit Hospitalization for Infectious Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
xaxis = list(title = "Date"),
yaxis = list(title = "Cumulative Deficit Hospitalization-Infectious",
tickformat = ","),
hoverlabel = list(bgcolor = "white",
font = list(family = "Arial", size = 12)),
hovermode = "closest")
direct_plotly#4
cumulative_deficit <- deficit_cumulative(fit_4x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
direct_plotly <- plot_ly(data = cumulative_deficit) %>%
add_ribbons(x = ~date,
ymin = ~observed - 2*sd,
ymax = ~observed + 2*sd,
name = "95% CI",
fillcolor = "rgba(180, 180, 180, 0.5)",
line = list(color = "transparent"),
hoverinfo = "none") %>%
add_lines(x = ~date,
y = ~observed,
name = "Observed",
line = list(color = "black", width = 2),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
add_markers(x = ~date,
y = ~observed,
name = "Data Points",
marker = list(color = "black", size = 4),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
layout(title = list(text = "Cumulative Deficit Hospitalization for infectious respiratory Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
xaxis = list(title = "Date"),
yaxis = list(title = "Cumulative Deficit Hospitalization-infectious respiratory",
tickformat = ","),
hoverlabel = list(bgcolor = "white",
font = list(family = "Arial", size = 12)),
hovermode = "closest")
direct_plotly#5
cumulative_deficit <- deficit_cumulative(fit_5x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
direct_plotly <- plot_ly(data = cumulative_deficit) %>%
add_ribbons(x = ~date,
ymin = ~observed - 2*sd,
ymax = ~observed + 2*sd,
name = "95% CI",
fillcolor = "rgba(180, 180, 180, 0.5)",
line = list(color = "transparent"),
hoverinfo = "none") %>%
add_lines(x = ~date,
y = ~observed,
name = "Observed",
line = list(color = "black", width = 2),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
add_markers(x = ~date,
y = ~observed,
name = "Data Points",
marker = list(color = "black", size = 4),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
layout(title = list(text = "Cumulative Deficit Hospitalization for injury Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
xaxis = list(title = "Date"),
yaxis = list(title = "Cumulative Deficit Hospitalization-injury",
tickformat = ","),
hoverlabel = list(bgcolor = "white",
font = list(family = "Arial", size = 12)),
hovermode = "closest")
direct_plotly#6
cumulative_deficit <- deficit_cumulative(fit_6x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
direct_plotly <- plot_ly(data = cumulative_deficit) %>%
add_ribbons(x = ~date,
ymin = ~observed - 2*sd,
ymax = ~observed + 2*sd,
name = "95% CI",
fillcolor = "rgba(180, 180, 180, 0.5)",
line = list(color = "transparent"),
hoverinfo = "none") %>%
add_lines(x = ~date,
y = ~observed,
name = "Observed",
line = list(color = "black", width = 2),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
add_markers(x = ~date,
y = ~observed,
name = "Data Points",
marker = list(color = "black", size = 4),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
layout(title = list(text = "Cumulative Deficit Hospitalization for pneumonia Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
xaxis = list(title = "Date"),
yaxis = list(title = "Cumulative Deficit Hospitalization-pneumonia",
tickformat = ","),
hoverlabel = list(bgcolor = "white",
font = list(family = "Arial", size = 12)),
hovermode = "closest")
direct_plotly#7
cumulative_deficit <- deficit_cumulative(fit_7x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
direct_plotly <- plot_ly(data = cumulative_deficit) %>%
add_ribbons(x = ~date,
ymin = ~observed - 2*sd,
ymax = ~observed + 2*sd,
name = "95% CI",
fillcolor = "rgba(180, 180, 180, 0.5)",
line = list(color = "transparent"),
hoverinfo = "none") %>%
add_lines(x = ~date,
y = ~observed,
name = "Observed",
line = list(color = "black", width = 2),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
add_markers(x = ~date,
y = ~observed,
name = "Data Points",
marker = list(color = "black", size = 4),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
layout(title = list(text = "Cumulative Deficit Hospitalization for renal Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
xaxis = list(title = "Date"),
yaxis = list(title = "Cumulative Deficit Hospitalization-renal",
tickformat = ","),
hoverlabel = list(bgcolor = "white",
font = list(family = "Arial", size = 12)),
hovermode = "closest")
direct_plotly#8
cumulative_deficit <- deficit_cumulative(fit_8x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
direct_plotly <- plot_ly(data = cumulative_deficit) %>%
add_ribbons(x = ~date,
ymin = ~observed - 2*sd,
ymax = ~observed + 2*sd,
name = "95% CI",
fillcolor = "rgba(180, 180, 180, 0.5)",
line = list(color = "transparent"),
hoverinfo = "none") %>%
add_lines(x = ~date,
y = ~observed,
name = "Observed",
line = list(color = "black", width = 2),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
add_markers(x = ~date,
y = ~observed,
name = "Data Points",
marker = list(color = "black", size = 4),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
layout(title = list(text = "Cumulative Deficit Hospitalization for respiratory Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
xaxis = list(title = "Date"),
yaxis = list(title = "Cumulative Deficit Hospitalization-respiratory",
tickformat = ","),
hoverlabel = list(bgcolor = "white",
font = list(family = "Arial", size = 12)),
hovermode = "closest")
direct_plotly#9
cumulative_deficit <- deficit_cumulative(fit_9x,
start = make_date(2020, 03, 01),
end = make_date(2024, 12, 28))
direct_plotly <- plot_ly(data = cumulative_deficit) %>%
add_ribbons(x = ~date,
ymin = ~observed - 2*sd,
ymax = ~observed + 2*sd,
name = "95% CI",
fillcolor = "rgba(180, 180, 180, 0.5)",
line = list(color = "transparent"),
hoverinfo = "none") %>%
add_lines(x = ~date,
y = ~observed,
name = "Observed",
line = list(color = "black", width = 2),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
add_markers(x = ~date,
y = ~observed,
name = "Data Points",
marker = list(color = "black", size = 4),
hoverinfo = "text",
text = ~paste("Date:", date,
"<br>Cumulative Deficit:", round(observed, 1),
"<br>Lower CI:", round(observed - 2*sd, 1),
"<br>Upper CI:", round(observed + 2*sd, 1))) %>%
layout(title = list(text = "Cumulative Deficit Hospitalization for other Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
xaxis = list(title = "Date"),
yaxis = list(title = "Cumulative Deficit Hospitalization-other",
tickformat = ","),
hoverlabel = list(bgcolor = "white",
font = list(family = "Arial", size = 12)),
hovermode = "closest")
direct_plotlyggsave("cumulative_deficit_plot_1.png", plot = cumulative_deficit_plot_1, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_2.png", plot = cumulative_deficit_plot_2, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_3.png", plot = cumulative_deficit_plot_3, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_4.png", plot = cumulative_deficit_plot_4, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_5.png", plot = cumulative_deficit_plot_5, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_6.png", plot = cumulative_deficit_plot_6, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_7.png", plot = cumulative_deficit_plot_7, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_8.png", plot = cumulative_deficit_plot_8, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_9.png", plot = cumulative_deficit_plot_9, width = 8, height = 6, dpi = 300)ggsave("cumulative_deficit_plot_1.png", plot = cumulative_deficit_plot_1, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_2.png", plot = cumulative_deficit_plot_2, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_3.png", plot = cumulative_deficit_plot_3, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_4.png", plot = cumulative_deficit_plot_4, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_5.png", plot = cumulative_deficit_plot_5, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_6.png", plot = cumulative_deficit_plot_6, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_7.png", plot = cumulative_deficit_plot_7, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_8.png", plot = cumulative_deficit_plot_8, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_9.png", plot = cumulative_deficit_plot_9, width = 5, height = 3, dpi = 300)# new 'category' column added
cardiovascular <- fit_x$detected_intervals %>%
mutate(category = "cardiovascular")
digestive <- fit_2x$detected_intervals %>%
mutate(category = "digestive")
infectious <- fit_3x$detected_intervals %>%
mutate(category = "infectious")
infectious_respiratory <- fit_4x$detected_intervals %>%
mutate(category = "infectious respiratory")
injury <- fit_5x$detected_intervals %>%
mutate(category = "injury")
pneumonia <- fit_6x$detected_intervals %>%
mutate(category = "pneumonia")
renal <- fit_7x$detected_intervals %>%
mutate(category = "renal")
respiratory <- fit_8x$detected_intervals %>%
mutate(category = "respiratory")
other <- fit_9x$detected_intervals %>%
mutate(category = "other")
combined_intervals <- bind_rows(
cardiovascular,
digestive,
infectious,
infectious_respiratory,
injury,
pneumonia,
renal,
respiratory,
other
)
combined_intervals <- combined_intervals %>%
select(category, everything())
combined_intervals <- combined_intervals %>%
mutate(across(c(where(is.numeric)), ~round(., 1)))
kable_table <- combined_intervals %>%
kable(
format = "html",
caption = "Detected Intervals by Disease Category",
align = c("l", rep("c", 9)),
col.names = c(
"Category", "Start Date", "End Date", "Observed Count Rate",
"Expected Count Rate", "SD Count Rate", "Observed Counts",
"Expected Counts", "Deficit", "SD"
),
row.names = FALSE
)
kable_table| Category | Start Date | End Date | Observed Count Rate | Expected Count Rate | SD Count Rate | Observed Counts | Expected Counts | Deficit | SD |
|---|---|---|---|---|---|---|---|---|---|
| cardiovascular | 2019-09-23 | 2019-12-02 | 8.1 | 8.7 | 0.1 | 12062 | 12840.1 | -778.1 | 113.3 |
| cardiovascular | 2020-01-20 | 2020-02-03 | 8.3 | 8.7 | 0.1 | 3343 | 3515.5 | -172.5 | 59.3 |
| cardiovascular | 2020-03-02 | 2021-02-15 | 7.3 | 8.9 | 0.0 | 49893 | 60938.8 | -11045.8 | 246.9 |
| cardiovascular | 2021-12-27 | 2022-01-31 | 8.3 | 9.0 | 0.1 | 6709 | 7290.1 | -581.1 | 85.4 |
| cardiovascular | 2023-05-01 | 2023-05-01 | 9.7 | 9.6 | 0.3 | 1303 | 1289.0 | 14.0 | 35.9 |
| digestive | 2019-09-15 | 2019-10-13 | 4.0 | 4.3 | 0.1 | 2684 | 2869.1 | -185.1 | 53.6 |
| digestive | 2020-03-01 | 2020-09-06 | 3.5 | 4.4 | 0.0 | 13267 | 16748.9 | -3481.9 | 129.4 |
| digestive | 2020-10-11 | 2021-03-14 | 3.7 | 4.2 | 0.0 | 11504 | 12987.3 | -1483.3 | 114.0 |
| digestive | 2021-12-05 | 2022-01-30 | 3.8 | 4.2 | 0.1 | 4556 | 5073.5 | -517.5 | 71.2 |
| digestive | 2022-05-01 | 2022-05-22 | 4.4 | 4.6 | 0.1 | 2353 | 2502.5 | -149.5 | 50.0 |
| digestive | 2022-09-18 | 2022-10-02 | 4.3 | 4.5 | 0.1 | 1722 | 1835.4 | -113.4 | 42.8 |
| infectious | 2020-05-03 | 2020-08-02 | 2.4 | 2.7 | 0.0 | 4495 | 5124.4 | -629.4 | 71.6 |
| infectious | 2020-10-04 | 2021-05-30 | 2.2 | 2.7 | 0.0 | 10517 | 12848.3 | -2331.3 | 113.4 |
| infectious | 2022-01-23 | 2022-02-27 | 2.5 | 2.8 | 0.1 | 2044 | 2250.7 | -206.7 | 47.4 |
| infectious | 2022-12-25 | 2023-01-22 | 2.6 | 2.8 | 0.1 | 1727 | 1898.7 | -171.7 | 43.6 |
| infectious | 2023-06-18 | 2023-07-02 | 2.5 | 2.8 | 0.1 | 1027 | 1119.5 | -92.5 | 33.5 |
| infectious | 2024-10-20 | 2024-12-15 | 2.5 | 2.8 | 0.0 | 3060 | 3354.2 | -294.2 | 57.9 |
| infectious respiratory | 2020-03-29 | 2021-05-30 | 0.1 | 0.9 | 0.0 | 1036 | 7153.0 | -6117.0 | 84.6 |
| infectious respiratory | 2021-10-31 | 2022-03-20 | 0.5 | 1.3 | 0.0 | 1333 | 3752.6 | -2419.6 | 61.3 |
| infectious respiratory | 2023-01-15 | 2023-02-26 | 0.7 | 1.3 | 0.0 | 616 | 1184.6 | -568.6 | 34.4 |
| infectious respiratory | 2023-04-09 | 2023-04-23 | 0.4 | 0.7 | 0.0 | 170 | 283.3 | -113.3 | 16.8 |
| infectious respiratory | 2024-11-03 | 2024-11-10 | 0.6 | 1.0 | 0.1 | 173 | 270.3 | -97.3 | 16.4 |
| injury | 2020-03-08 | 2020-09-20 | 6.3 | 7.8 | 0.0 | 24783 | 30324.6 | -5541.6 | 174.1 |
| injury | 2020-10-04 | 2021-02-14 | 6.6 | 7.5 | 0.1 | 17808 | 20166.2 | -2358.2 | 142.0 |
| injury | 2021-03-07 | 2021-04-11 | 7.1 | 7.6 | 0.1 | 5752 | 6123.8 | -371.8 | 78.3 |
| injury | 2021-05-09 | 2021-05-16 | 7.7 | 8.0 | 0.2 | 2082 | 2147.9 | -65.9 | 46.3 |
| injury | 2021-07-25 | 2021-08-08 | 8.0 | 8.4 | 0.1 | 3223 | 3381.7 | -158.7 | 58.2 |
| injury | 2021-12-19 | 2022-01-16 | 7.1 | 7.7 | 0.1 | 4807 | 5171.7 | -364.7 | 71.9 |
| injury | 2022-07-31 | 2022-08-21 | 8.3 | 8.7 | 0.1 | 4453 | 4668.8 | -215.8 | 68.3 |
| pneumonia | 2020-04-26 | 2021-07-04 | 1.5 | 2.6 | 0.0 | 12636 | 22211.3 | -9575.3 | 149.0 |
| pneumonia | 2021-09-05 | 2022-05-22 | 1.9 | 2.8 | 0.0 | 9964 | 14396.8 | -4432.8 | 120.0 |
| pneumonia | 2022-07-17 | 2022-08-28 | 1.7 | 2.0 | 0.0 | 1597 | 1931.8 | -334.8 | 44.0 |
| pneumonia | 2022-10-30 | 2022-11-13 | 2.6 | 2.9 | 0.1 | 1035 | 1163.3 | -128.3 | 34.1 |
| pneumonia | 2023-01-22 | 2023-02-19 | 2.4 | 3.0 | 0.1 | 1621 | 2005.2 | -384.2 | 44.8 |
| pneumonia | 2023-04-23 | 2023-05-21 | 2.3 | 2.7 | 0.1 | 1545 | 1792.5 | -247.5 | 42.3 |
| pneumonia | 2024-01-28 | 2024-02-11 | 2.6 | 3.0 | 0.1 | 1058 | 1201.7 | -143.7 | 34.7 |
| renal | 2020-03-08 | 2020-07-12 | 1.8 | 2.3 | 0.0 | 4551 | 5929.0 | -1378.0 | 77.0 |
| renal | 2020-09-06 | 2020-09-27 | 2.2 | 2.3 | 0.1 | 1182 | 1260.9 | -78.9 | 35.5 |
| renal | 2020-11-01 | 2020-12-20 | 1.9 | 2.2 | 0.0 | 2074 | 2355.1 | -281.1 | 48.5 |
| renal | 2021-01-03 | 2021-02-28 | 2.0 | 2.2 | 0.0 | 2382 | 2723.9 | -341.9 | 52.2 |
| renal | 2022-01-02 | 2022-01-02 | 2.1 | 2.3 | 0.1 | 286 | 314.5 | -28.5 | 17.7 |
| renal | 2022-05-15 | 2022-05-29 | 2.4 | 2.6 | 0.1 | 955 | 1053.3 | -98.3 | 32.5 |
| respiratory | 2020-03-08 | 2021-05-30 | 2.7 | 4.2 | 0.0 | 23417 | 36478.4 | -13061.4 | 191.0 |
| respiratory | 2021-11-21 | 2022-03-27 | 3.6 | 4.4 | 0.0 | 9293 | 11363.2 | -2070.2 | 106.6 |
| respiratory | 2023-04-23 | 2023-05-07 | 4.2 | 4.7 | 0.1 | 1691 | 1897.1 | -206.1 | 43.6 |
| other | 2020-03-08 | 2021-05-23 | 35.6 | 41.6 | 0.1 | 307080 | 358695.6 | -51615.6 | 598.9 |
| other | 2021-08-22 | 2021-10-10 | 40.6 | 42.8 | 0.2 | 43766 | 46098.9 | -2332.9 | 214.7 |
| other | 2021-11-21 | 2022-02-13 | 36.0 | 40.3 | 0.2 | 63018 | 70484.7 | -7466.7 | 265.5 |
| other | 2022-03-06 | 2022-04-10 | 41.3 | 43.1 | 0.2 | 33364 | 34828.8 | -1464.8 | 186.6 |
| other | 2022-04-24 | 2022-05-08 | 42.6 | 44.3 | 0.3 | 17206 | 17908.5 | -702.5 | 133.8 |
library(kableExtra)##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rowslibrary(webshot2)
# Create styled HTML table
html_table <- combined_intervals %>%
kable("html",
escape = FALSE,
caption = "Detected Intervals by Disease Category",
align = c("l", rep("c", ncol(combined_intervals) - 1)),
col.names = c(
"Category", "Start Date", "End Date", "Observed Count Rate",
"Expected Count Rate", "SD Count Rate", "Observed Counts",
"Expected Counts", "Deficit", "SD"
)) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE,
position = "center")
# Save as HTML and then convert to PNG
save_kable(html_table, "combined_intervals_table.html")
webshot("combined_intervals_table.html", "combined_intervals_table.png", vwidth = 1200, vheight = 800)